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ABSTRACT 



Narrow-band Lya line and broad-band continuum have played important roles in 
the discovery of high-rcdshift galaxies in recent years. Hence, it is crucial to study the 
Q ■ radiative transfer of both Lya and continuum photons in the context of galaxy for- 

mation and evolution in order to understand the nature of distant galaxies. Here, we 
C/3 . present a three-dimensional Monte Carlo radiative transfer code, All-wavelength Ra- 

^ ■ diative Transfer with Adaptive Refinement Tree (ART^), which couples Lya line and 

multi-wavelength continuum, for the study of panchromatic properties of galaxies and 
■ interstellar medium. This code is based on the original version of Li et al., and features 

^ ' three essential modules: continuum emission from X-ray to radio, Lya emission from 

both recombination and coUisional excitation, and ionization of neutral hydrogen. The 
0^ ■ coupling of these three modules, together with an adaptive refinement grid, enables a 

OO I self-consistent and accurate calculation of the Lya properties, which depend strongly 

. on the UV continuum, ionization structure, and dust content of the object. More- 

Q^ ■ over, it efficiently produces multi-wavelength properties, such as the spectral energy 

I distribution and images, for direct comparison with multi-band observations. 

As an example, we apply ART^ to a cosmological simulation that includes both 
star formation and black hole growth, and study in detail a sample of massive galaxies 
at redshifts z = 3.1 — 10.2. We find that these galaxies are Lya emitters (LAEs), 
whose Lya emission traces the dense gas region, and that their Lya lines show a 
' shape characteristic of gas inflow. Furthermore, the Lya properties, including photon 

M ■ escape fraction, emergent luminosity, and equivalent width, change with time and 

environment. Our results suggest that LAEs evolve with redshift, and that early LAEs 
such as the most distant one detected at z ~ 8.6 may be dwarf galaxies with a high 
star formation rate fueled by infall of cold gas, and a low Lya escape fraction. 

Key words: radiative transfer ~ line: profiles ~ ISM: dust, extinction - galaxies: 
evolution - galaxies: formation - galaxies: high-redshift 



1 INTRODUCTION 

The hydrogen Lya emission is one of the strongest lines 
in the UV band. It provides a unique and powerful 
tool to search for d i stant galaxies, as first suggested by 
[Partridge fc PeeblesI (|l967l '). and evidenced by the suc- 
cessful detection of hundreds of Lya emitting galaxies, 
or Lya emitters (L AEs), at redshif t s z > 3 over the 
past decade (e.g., IHu fc McMahonI 1 19961 : ICowie fc Hul 
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19981: iHu. Cowie. fc McMahonI Il998l: ISteidel et all 200] 



Rhoads et all l2000l: IFvnbo. MoUer fc ThomsenI 
Hu et all 12*00^; iFvnbo et all l2003l: iRhoads ct al.' 



Ouchi et al.] 120031: iKodaira et all 12003 
■ 12004'; Dawson et al 



Hu et al 



2004 



20011 : 

200j 

Maicr ct al. 200^ 
Malhotra fc Rhoadd 



20051: iKashikawa et all l2006l: IShimasaku et all 


20061 


Ive et al. 20061: Hu & Cowid 


20061: Gronwall et all 


200d 


Cubv et all 


20071: IStark ct al 


l2007l: .Nilsson et al.[ 


20071 


Ouchi et al.l 


20081: IWillis. Courbin. Kneib. fc Minnitil 


2008| 



Ota et al. 20081: IHu et all 



2010l : 



Ouchi et all 



2010l ). 
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recently, iLehnert et al.l (|201Cf ) have discovered the most 
distant galaxy at z = 8.6 using the Lya line. 

To obtain a full picture of the physical properties 
of the high-redshift LAEs, efforts have been made re- 
cently to_observe_theseobiect^ c ontinuum 
(e.g., iGawiser et al ] I2OO6I: iGronwall et all l2007l: I Lai et al.l 



2007i:lNilsson et al.1 2OT7I: Pirzkal et al.ll2C)07l:lLai et al.||20o' 



Ouchi et al.'2008': Pcntcricci ct al, 2009; Ono ct al. '20f0a|L 
Haves et al..,20fQ : .Finkclstcin ct al. 2011 ; Nilsson fc Molleil 
20 111 ). One remarkable result from these multi-wavelength 



surveys is that, LAEs appear to evolve with cosmic time: 
the LAEs at high redshifts z > 5 tend to be smaller and 
younger galaxies, and have less dust extinction and higher 
equivalent widths (EWs) than their counterparts at lower 
redshifts z ~ 3. 

Despite the rapid progress in the Lya detections, 
the origin and nature of these distant LAEs, however, 
remain largely unknown, because the radiative transfer 
(RT) of Lya line is a complicated process. It depends 
on the geometry, kinematics, ionization state, and dust 
content of the surrounding medium, as well as the stellar 
population and other ionizing sources of the galaxy. To 
date, there have been a number of theoretical studies on the 
resonant scatter ing of Lya pho t ons, a nalyticall y in simple 



geom e try ( e. g., HummeJ Il962l; Adams 197^; iHarringtonl 



Il973l . Il974i : iNeufeldl Il99d. Il99ll: iLoeb fc Rvbickil 



and numerically (e.g . T^Auej" Il968l: lAverv fc House! Iigei; 



Ahn. Lee, fc Ledl200oi.[200ll.l2002l:IZheng fc Miralda-Escudel 



Diikstra. Haiman. fc Spaanj I2OO6I: iHansen fc Oh 
Tasitsiomil I2OO6I: I mme. Schaerer. fc Masellil 
Laursen fc Sommer-LarsenI l2007l : 

Laursen. Razoumov. fc Sommer-LarsenI l2009al: 

Pierleoni. Maselli. fc Ciardil l2009l: iFaucher-Giguere et al.l 



20ld : IZheng et al.ll201(]| . I2OIII : ISchaerer et al.ll201lh . Among 
the numerical codes, Monte Carlo method is widely 
employed to solve the Lya RT, owing to the irregular 
geometry and inhomogeneous distribution of the interstellar 
medium (ISM), and the slow convergence of the numerical 
techniques at high optical depths. By applying Lya RT 
codes to cosmological simulations, several aspects of the 
LAEs have been investigated by various groups, including 
obser ved properties of LAEs at z ~ 5.7 (jZheng et al.ll201o[ 
I2OIII). escape fraction of Lya photons from hig h-z galaxies 
l|Laursen. Sommer-Larsen. fc AndersenI l2009bl ). effects of 



ioniz ation on the emergent Lya spectra ( Pierleoni et al.l 
I2OO9I ) and mechanism o f high EWs in dusty multi-phase 
ISM (|Hansen fc Ohll2006l '). However, most of these previous 
works have focused only on the Lya line, and were therefore 
limited to some aspect of the Lya properties, such as flux 
and line profile. 

On the other hand, there has been an ar- 
ray of RT codes in continuum over the last sev- 
eral decades, u si ng ei the r ray-tracing methods (e.g.. 



Rowan- RobinsonI 1980l: lEfstathiou fc Rowan-RobinsonI 



I990I: [Steinacker et all l2003l: iFolini et all l2003l : 



Steinacker. Bacmann. fc Hennina 2009), or Monte Carlo 



Lefevre. Daniel. & Bcrgcat ,1983.: 


.Whitncv & Hartmann 


19921: Witt, Thronson. & Caouanol 


1992: Code fc Whitney 


I995I: iLooez, Mekarnia, & Lefevre 


I1995I: ILucvI Il99d: 


Wolf. Henning. fc Stccklum 19991: Bianchi. Davies. fc Alton 


2OOOI: Biorkman & WoodI I2001I: 


Whitnev et all 20031: 



Harries et al. 


Li et al. 


2008 


Li et al. 


(200 



20041: [jonssoni I2OO6I: IPirite et all I2OO6I : 



developed a multi-wavelength RT code that 
employs radiative equilibrium and an adaptive refinement 
grid, which reproduced the observed dust properties of the 
most distant quasars detected at z ^ 6 when applie d to the 
hydrodynamic quasar simulations of iLi et al] (|2007l ). 

However, most of these developments focused mainly on 
radiation transport in dusty environments, in particular the 
absorption and re-emission by dust, of local star-forming re- 
gions. In order to explain the diverse properties of LAEs and 
their evolution with redshift, it is critical to couple Lya with 
continuum, in the context of galaxy formation and evolution. 
Several aspects of such an appr oach have been attempted in 
previous studies. For examples. iPierleoni et all (|2009l ) incor- 
porated prfr;Coni£uted_grid of Lya RT into ionization calcu- 
lations, iTasitsiomil (|2006l) considered a combination of Lya 



and ion ization RT, Laursen. Sommer-Larsen. fc AndersenI 



(|2009bl l studied the effect of ionization and dust on transfer 
of Lya photons, while ISchaerer et al.l (|201lf ) included dust 
and non-ionizing UV continuum in the Lya RT calculations. 
An ultimate goal would be to tackle the ionization of hy- 
drogen, non- ionizing continuum, interstellar dust, and Lya 
propagation and scattering simultaneously. 

Furthermore, the implementation of Lya emission 
should include two major generation mechanisms, the re- 
combination of ionizing photons, and coUisional excitation 
of hydrogen gas. In addition, relevant ionizing sources such 
as stars, active galactic nucleus (AGN), and UV background 
radiation should be included. Last, but not the least, the RT 
code should incorporate an adaptive refinement grid in or- 
der to cover a large dynamical range and resolve dense gas 
regions where star for mation takes place an d prominent Lya 
emission comes from (|Laursen et al.l^2009a^ . 

In this work, we present a new, 3-D Monte Carlo RT 
code. All-wavelength Radiative Transfer with Adaptive Re- 
finement Tree (ART^), which couples Lya line and multi- 
wavelength continuum, for the study of panchromatic prop- 
erties of galaxies and ISM. (Note that the "coupling" in 
this paper means the RT of Lya with other continuum pho- 
tons, not the coupling between RT and hydrodynamics.) Our 
code improves over the original continuum-only version of 
iLi et al.l (|2008l ), and features three essential modules: contin- 
uum emission from X-ray to radio, Lya emission from both 
recombination and coUisional excitation, and ionization of 
neutral hydrogen. The coupling of these three modules en- 
ables a self-consistent and accurate calculation of the Lya 
and multi-wavelength properties of galaxies, as the equiva- 
lent width of the Lya line depends on the UV continuum, 
and the escape fraction of Lya photons strongly depends on 
the ionization structure and the dust content of the object. 
Moreover, the adaptive refinement grid handles arbitrary 
geometry and efficiently traces inhomogeneous density dis- 
tribution in galaxies and ISM. Furthermore, this code takes 
into account radiation from both stars and accreting black 
holes, so it can be used to study both galaxies and AGNs. 
ART^ can produce a number of observables, including spec- 
tral energy distribution (SED) from X-ray to radio, multi- 
band fluxes and images, Lya emission line and its EW, which 
can be directly compared to real observations. It has a wide 
range of applications, and can be easily applied to simula- 
tions using either grid-based or smoothed particle hydrody- 
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namic (SPH) codes by converting the simulation snapshot 
to the grid structure of ART^. 

The paper is organized as follows. In §2, we describe the 
ART^ code, which includes implementation of three mod- 
ules of Continuum, Lya line, and Ionization of hydrogen, 
the adaptive refinement grid, and the dust model. In §3, 
we present the application of ART^ to a SPH cosmological 
simulation, and study in details the Lya emission and multi- 
band properties from individual galaxies at redshifts z=8.5, 
6.2, and 3.1 from the simulation. We discuss in §4 contri- 
bution to Lya emission from AGNs, stars, and excitation 
cooling from gas accretion in our model, and effects of nu- 
merical resolutions on the RT calculations, and summarize 
in 55. 



2 ART^: ALL- WAVELENGTH RADIATIVE 
TRANSFER WITH ADAPTIVE 
REFINEMENT TREE 

ART'^ is a 3-D, Mo n te Ca rlo RT code based on the original 
version of iLi et al.l (|2008h , which included the continuum 
emission from X-ray to radio, and an adaptive refinement 
grid. We have added two more modules, Lya emission from 
both recombination and coUisional excitation, and ioniza- 
tion of neutral hydrogen, to the current version and couple 
them self- consistently. Th e continuum part was described in 
details in iLi et all (|2008l 'l. Here we briefly outline the Con- 
tinuum procedure, and focus on the new implementations of 
Lya line and Ionization. 



2.1 Continuum Radiative Transfer 



The C ontinuum module in ART^ was developed in lLi et al] 
(|2008|). which adopted the r adiative equilibrium algorithm 
of iBiorkman fc WoodI (|200ll '). In dusty environments, ab- 
sorbed radiation energy by dust is re-emitted as thermal 
emission. The re-emitted spectrum depends on the tem- 
perature of the dust, which is assumed to be in thermal 
equilibrium with the radiation field. The radiative equilib- 
rium is ensured by performing the Monte Carlo transfer it- 
eratively until the dust temperature distribution converges, 
which can be computationally expensive. To accelerate the 
calculation, ART^ uses th e "immediate reemission" scheme 
jBiorkman fc WoodI [iooH ) . in which the dust temperature 
is immediately updated on absorption of a photon packet, 
and the frequency of re-emitted photons are sampled from a 
spectrum that takes into account the modified temperature. 
The temperature is determined by the balance between the 
stacked energy absorbed by dust in the cell and the thermal 
emission from them. The emitted energy in a cell in the time 
interval At is 



Ef 



in At J dVi J pK^B^{T)di' 
4-KAtnp{T^)B{T,)mi , 



(1) 



where Kp = J n^B^ di// B is the Planck mean opacity, B = 
aT^ j-K is the frequency integrated Planck function, and rrii 
is the dust mass in the cell, the subscript i indicates the i'th 
cell. 

Solving the balance between the absorbed and emitted 



energy, we obtain the dust temperature as follows after ab- 
sorbing Ni packets. 



'iN^K,p(Ti)mi 



(2) 



where A''.^ is the total number of photon packets in the sim- 
ulation, and L is the total source luminosity. Note that 
because the dust opacity is temperature-independent, the 
product Kp{Ti)aTi increases monotonically with tempera- 
ture. Consequently, Ti always increases when the cell absorbs 
an additional packet. 

The added energy to be radiated owing to the tempera- 
ture increase AT is determined by a temperature-corrected 
emissivity Aj^ in the following approximation when the 
temperature increase, AT, is small: 



Aju ~ KvpAT 



dB^T) 
dT 



(3) 



The re-emitted packets, which comprise the diffuse ra- 
diation field, then continue to be scattered, absorbed, and 
re-emitted until they finally escape from the system. This 
method conserves the total energy exactly, and does not re- 
quire any iteration as the emergent SED, uL,, — k,^B^{T), 
corresponds to the equilibrium temperature distribution. 



2.2 Lya Line Transfer 

Hydrogen Lya photon corresponds to the transition between 
the n = 2 and n — 1 levels of a hydrogen atom. It is 
the strongest Hi transition. The RT of Lya photons is de- 
termined by Lya resonant scattering, dust absorption and 
scattering, and ionization state of the medium. The pro- 
cess is highly complicated in galaxies owing to the com- 
plex geometry and gas distribution. In order to acceler- 
ate the numerical convergence of the RT process, Monte 
Carlo metho d has been commonly used in a number of Lya 



codes (e.g., [Zh cng fc Miralda-Escudi |2002|: 



200d:lTasits"iomi2006:,yerhainme et al.ll2006l: 



Diikstra et al 



Laursen et al 



2009al : iPierleoni et al.l |2009| : iFaucher-Giguere et al.l l201(il ) 

Our implementation of Lya line transfer adopts the Monte 
Carlo method, and the major improvements over many of 
these codes are that, it is coupled with ionization and multi- 
wavelength continuum which enables a self-consistent and 
accurate calculation of the Lya properties, and is incorpo- 
rated with a 3-D adaptive-mesh refinement grid which effi- 
ciently handles arbitrary geometry and inhomogeneous den- 
sity distribution. Moreover, we treat both recombination and 
coUisional excitation for Lya emission. 



2.2.1 Propagation and Scattering of Lya Photons 

The optical depth Tv{s) of a Lya photon with frequency u 
traveling a path of length s is determined by 



r„(s) = 



?l(Vj|)cr„ dV\\dl, 



(4) 



where ) is the number density of neutral hydrogen gas 
with parallel velocity component Vn , Oi, is the scattering 
cross section as a function of frequency. In the rest frame of 
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the hydrogen atom, Oi, takes the form 



Ov — fi: 



Ai^L/27r 



eC {U - H,)'' + (AfL/2)2 



(5) 



where /12 — 0.4162 is the Lya oscillator strength, vq = 
2.466 X 10^^ Hz is the line-center frequency, Aul = 9.936 x 
10^ Hz is the natural line width, and the other symbols have 
their usual meaning. Assuming a Maxwellian distribution for 
the thermal velocity of the encountering atoms, the resulting 
average cross section is 



where 



/i: 



H{a,x) = — 



H{a,x), 



dy 



(6) 



(7) 



is the Voigt function, Az/d = [2kBT/{TnpC'^)]^^^uo is the 
Doppler width, x — {1/ ~ uo)/Ai'o is the relative frequency 
of the incident photon in the laboratory frame, and a = 
Ai/L/(2At'D) is the relative line width. 

Some previous work approximated the Voigt function 
with a Gaussian fitting in the core and a power law fitting 
in the wing. However, since this approximation causes a rel- 
ative error of a few te ns per cent at the transit region, we 
use the analytical fit of iTasitsiomil ()2006l ). 



H{a, x) 



■ q^TT + e 



where 



<? = 



for C < 

{^ + ^)vuhT)Tl(0 forC>0, 



(8) 



(9) 



with C = {x^ - 0.855)/(a;2 + 3.42) and U(Q = 5.674^* - 
9.207C^ -h4.421C^ -^0.1117C. 

This approximation fits the Voigt function well for all 
frequencies, and the relative error is always less than 1 per- 
cent above a temperature of 2K. 

The optical depth of dust absorption and scattering is 
estimated as. 



dr — nmo'{i')ds + mdced{i')ds 



(10) 



where rud is the dust mass and ad = Qd.abs + Qd.sca is the 
mass opacity coefficient of absorption and scattering. When 
the stacked optical depth through passing cells achieves an 
optical depth determined by Equation (|10p . the encounter- 
ing medium is chosen using a random number by comparing 
with the fraction, nHi<j{i')/ (nHio-(i/) 4- mdctd{i^))- 

When Lya photons are scattered by neutral hydrogen 
atoms, the frequency in laboratory frame is changed depend- 
ing on the velocity components of the atoms and the direc- 
tion of the incidence and the scattering. Then the velocity 
components of the directions perpendicular to the incident 
direction will follow a Gaussian distributi on, and can be gen - 
erated by a simple Box-MuUer method (|Press et al.lll993 ). 
However, the parallel component of the velocity depends on 
relative frequency x of the incident photon owing to the res- 
onance nature of the scattering. The probability distribution 
of the parallel component is drawn from 



/hi) = 



TTH(a, x) (x — it||)2 + 



(11) 



where uy = V||/wth is velocity of the parallel component 
normalized by thermal velocity Wth (hereafter u means nor- 
malized veloci ty by Dth). To follow this distrib ution, we use 
the method of IZheng fc Miralda-Escudel (|2002l '). 

When a photon has a frequency x < Xcw in the 
optically-thick cell, where Xcw is the boundary between the 
core and the wing of the Voigt profile, i.e., where /^/n = 
a/nx'^, it cannot travel long distance and is confined in the 
cell with numerous scattering. Only when it has a large x 
from scattering by high velocity atoms can it escape from the 
cell and travel long distance. The photon usually experiences 
many scatterings (A'^ > 10^) before it moves into the wing, it 
is therefore extremely computation costly to trace this pro- 
cess for all photon packets. In order to speed up the calcu- 
lation, in particular to avoid the huge number of scatterings 
in the cor e , we use a core skipping method developed by 
lAhn et al] (I2OO2I') and follow the pr ocedure of lDiikstra et al.l 
(|2006l ) and iLaursen et al.l (1200931 ). It artificially push the 
photon in the wing by scattering with atoms which have 
high velocity in the direction perpendicular to the incident 
direction. The velocity of the perpendicular component u± 
is generated from a truncated Gaussian, 



wx,i = (a;crit ~ InRi) ' cos27rR2 
u±,2 = {x%i^ — lnRi)^''^sin27rR2. 



(12) 



where Ri and R2 are two univa riates, and we use the criti- 
cal frequency Xait introduced in lLaursen et al.l (|2009al ), i.e., 
a;crit = 0.02e«('"''^o''' where (C,x) = (0.6, 1.2) for am 60, 
or (Ci x) = (1.4,0.6) for aro > 60. This acceleration scheme 
can reduce the calculation time by several orders of magni- 
tude, and can produce a line profile which agrees well with 
analytical solutions in a static slab. 

The frequency after scattering depends on the direction 
in which the photon is scattered. The direction is given by 
the phase function, W{9) oc 1 + ^cos'^6', where 9 is the angle 
between the incident direction rii and the outgoing direction 
rif, and R/Q is the degree of polarization. The ratio R/Q 
becomes 3/7 for the 2P3/2 state in 2; < a^ rw (iHamilton [T94Cll ) ■ 
and 1 for the scattering in the wing jStenflcil 198C ) . The 
transition from 2P1/2 in the core, together with the core 
skipping scheme, results in isotropic scattering. 

The difference in the phase function does not affect the 
LyQ properties such as the escape fraction /esc, the emer- 
gent line profile and luminosity, even if a single phase func - 
tion is used for all scatterings (e.g., iLaursen et al.ll2009al l. 
Therefore, we simply assume an isotropic scattering for all 
the scattering process in this work. We should point out, 
however, that the phase function b ecomes important in th e 
polarization effect of Lyo? line (e.g.. iDiikstra fc LoebllioOSl l. 
and so care must be taken in dealing with that. 

In the scattering process, the final frequency in the lab- 
oratory frame is then 



■ U|l + n/ • u — g(l — Yii ■ lif) 



(13) 



where the factor g — h pfo/mHCVth takes into account t he re- 
coil effect (Field 1959|; IZheng fc Miralda-EscudelbOO^ ) with 
hp being the Planck constant. We trace the Lya photon 
packet until it is absorbed by dust or escapes from a galaxy 
(e.g., when the photons move out of the calculation box typ- 
ically 10 times larger than the virial radius of a galaxy). 
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2.2.2 Lya Emissivity 

Lya emission is generated by two major mechanisms: re- 
combination of ionizing photons and collisional excitation 
of hydrogen gas. 

• Recombination: Ionizing radiation from stars, AGNs 
and UVB can ionize the hydrogen gas in galaxies. The colh- 
sion by high-temperature gas can also ionize the hydrogen. 
The ionized hydrogen atoms then recombine and create Lya 
photons via the state transition 2P — !> IS. The Lya emissiv- 
ity by the recombination is 



rcc J. 1 

= JaOLBILVancnmi, 



(14) 



where qb is the case B recombination coefficient, and 
fa is the average number of Lya photons produced 
per case B recombi nation. Here we use as derived in 
iHui fc Gnedliil (| 19971 ). Since the temperature dependence 
of fa is not strong, fg = 0.68 is assumed everywhere 
lOsterbrock fc Ferlandll2006l ). The product hva is the en- 
ergy of a Lya photon, 10.2 eV. 

• Collisional Excitation: High temperature electrons 
can excite the quantum state of hydrogen gas by the colli- 
sion. Due to the large Einstein A coefficient, the hydrogen 
gas can occur de-excitation with the Lya emission. The Lya 
emissivity by the collisional excitation is estimated by 



coll ^ 
^a = (^LyctWenHI, 



(15) 



where 



Chyc 

3.7 



is the collisional excitation coefficient. 



X 10" 



(|Osterbrock fc FerlandllioO^ ') 



pi-hUa/kT)T 



-1/2 



ergs s 



Once the ionization structure have been determined (see 
Section r2.3[l . we estimate the intrinsic Lya emissivity in each 
cell by the sum of above Lya emissivity, ea = eU*^ + £q°''- 
The excitation cooling dominates at Tgas ~ 10*~* K, but 
becomes s maller than the recombination cooling at Tgas > 
10*^ K fe.g.. lFaucher-Giguere eralll2010l ). 

There is a large ambiguity in the estimation of the ex- 
citati on cooling from multi-compo nent ISM in SPH simula- 
tions iFaucher-Giguere et al.ll201Cll ). However, in the present 
work, we find that the power from stars and AGNs is al- 
ways larger than or comparable with the excitation cooling. 
Hence, even the largest possible cooling rate is still sub- 
dominant to the nebular Lya emission (i.e., coming from 
Hll regions around young hot stars). In this paper, we esti- 
mate the ionization rate and Lya emissivity for each cell by 
using the mixed physical quantities of multi-component ISM 
at first, then reduce it by weighting with the mass fraction 
of cold gas /cold, i.e., e^"" = eJi°o x (1 - /cold), where /coid 
can reach ~ 0.9 depending on the situation. 



2.2.3 Test Calculations 

To test our implementation of the Lya RT, we perform some 
standard tests against analytical solutions, as well as other 
numerical results in the literature. 

1. Neufield Test 

As the first test, we carry out the RT calculation in a 
dust-free slab of uniform gas. T he uniform slab was analyt- 
ically studied bv lNeufeldl (|l990l ) in the optically thick limit. 



0.01 



0.008 



0.006 



0.004 



0.002 - 




Figure 1. The Neufield Test of emergent LyQ profiles of 
monochromatic lino radiation emitted from a dust-free uniform 
slab. The black lines are our simul ation results, a nd the color 
lines represent analytical solutions bv lNeufeldl l|l990h for different 
optical depths: red — tq = 10*, blue - tq = 10^, and green — 
TO = 10^. The temperature of gas is 10 K. 



In this case, the emergent Lya line profile is given by 

, . Ve 1 

J(±T0,X) = 



24 



(16) 



where tq is the optical depth at the line center from mid- 
plane to the boundary of the slab and Xinj is the injection 
frequency x. 

Figure [T] shows our test result of emergent spectra from 
a dust-free gas slab with a temperature of T = 10 K and a 
central plane source, in comparison with analytical solutions 
for optical depths tq = 10^, 10^ and 10^, respectively. The 
width of the peaks becomes larger with larger to, in agree- 
ment with the dependence of the optical depth on frequency. 
Our simulation results agree very well with the analytical 
solutions, and the agreement becomes better with higher tq. 

The dust absorption is crucial in the study of es- 
cape fra c tion, flux, image and proflle of Lya from galaxies. 
iNeufeldl (| 19901 ) has derived an approximate expression for 
the escape fraction of Lya photons in a static dusty slab, 

/esc = r ^ (17) 



cosh 



,1/3, 



where Ta is the absorption optical depth of dust, and ^ = 
V^/C'T^''^^, with ^ ~ 0.525 being a fitting parameter. This 
solution is valid for extremely optically thick which arc ^ 
10'^, and in the limit of (aro)^^^ S> Ta. It suggests that the 
escape fraction decreases rapidly with increasing (aro)'^^^Ta. 

In figure [21 we compare our simulation results of Lya 
escape fraction in a dusty slab with the analytical solutions 
of Ncufcld (199^. As is shown, our results agree with the 
analytical curve very well. 

2. Loeb & Rybicki Test 
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Figure 2. The Neufield Tost of escape fraction of Lya photons 
in a dusty slab as a function of (a tq)^^^ Tabs- The solid lino is an 
approximate analytical solution of iNeufeldl lll990l ). and the open 
circles are our simulation results. 



Figure 3. The Loeb & Rybicki Test of mean intensity profiles of 
a monochromatic Lya source in a uniformly expanding medium 
at various frequenc i es. Th e dotted lines are analytical solutions 
by iLoeb fc Rvbickil l|l999l ) in the diffusion limit, and the color 
lines are our simulation results. 



iLoeb fc Rvbickil fTgoi) analytically derived the inten- 
sity field in a spherically symmetric, uniform, radially ex- 
panding neutral hydrogen cloud surrounding a central point 
source of Lya photons. No thermal motion is included 
{T = Q K). In the difi'usion limit, the mean intensity J(f , D) 
as a function of distance from the source f and frequency 
shift V is given by 

with v — {uo — v)lvi,^ where Vi, is the comoving frequency 
shift of Lyof at which the optical depth to infinity is unity, 
and f = r/r^ is the scaled distance, where is the physical 
distance at which the Doppler shift from the source due to 
the Hubble- like expansion equals Vi,. 

We use a s imulation setup similar to that of 
ISemelin. Combes, fc Back (2007i) for the test. The results 
are present ed in Figure [3] in c ompa rison with the analytical 
solutions of iLoeb fc Rvlaickil (| 19991 '). Our simulation results 
are very close to those fro m previous tests (e.g.. pKisitsiomil 
l2006l : ISemefin et al.|[2007l ). and are in good agreement with 
the analytical solutions in regimes where the diffusion limit 
is valid. However, the simulations diverge from the analytic 
solutions at larger f, where the assumption of optically thick 
is no longer valid. 

3. Bulk Motion Test 

The bulk motion of gas affects the escape of 
Lyg photons as it decreases the effecti ve optical depth 
(|Diikstra et al.|[2006l : iLaursen et "al]|2009al ). In practice, the 
bulk speed of the gas to the center can be ~ 10— 1000 km s^^ 
for infiow caused by gravitation, or outfiow by supernovae. 
The relative velocity can decrease the optical depth sig- 
nificantly. Unfortunately, there is no analytical solution of 



emer gent spectrum f or mov ing medium of T 7^ K. How- 
ever, |Laurse^^et^l] l|2009al ) calculated the emergent spec- 
tra of an isothermal and homogeneous sphere undergoing 
isotropic expansion or contracti on. Here, we follow the pro- 
cedure of iLaurseri et al.l (|2009al ) to set up the test. The ve- 
locity Vbuik(r) of a fluid element at a distance r from the 
center is set to be 

Vbuik(r)=Hr, (19) 

where the Hubble-like parameter "H is fixed such that the 
velocity increases linearly from in the center to a maximal 
absolute velocity Umax at the edge of the sphere {r — R): 

n = (20) 

with iimax positive for an expanding sphere, or negative for 
a collapsing one. 

Figure|4]shows the our test results of the emergent spec- 
trum from an isothermal and homogeneous sphere under- 
going isotropic expansion with different maximal velocities 
«max at the edge of the sphere. As Umax increases, the blue 
wing is suppressed while the red wing is broadened. This is 
because in an outfiowing sphere, the photons in the blue 
wing have a higher probability of encountering hydrogen 
atoms and being scattered than those in the red wind, which 
can escape more easily. This plot suggests that the peak of 
the profile is pushed further away from the line center for 
increasing Umax- However, above a certain threshold value, 
the peak moves back toward the center again owing to the 
decrease of optical depth caused by the fast bulk motion. 
Our result s show a good agreem ent with the previous simu- 
lations by iLaursen et al.l (|2009al '). 
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Figure 4. Bulk Motion Tost of emergent Lyo spectrum from an 
isothermal and homogeneous sphere undergoing isotropic expan- 
sion with different maximal velocities at the edge of the sphere. 
The initial conditions are the same as those used for Figure 8 in 
iLaursen et alj l l2009al ). with a temperature of the sphere T = 10* 
K, and the hydrogen column density TVui = 2 X 10^^ cm~^. 



2.3 Ionization Radiative Transfer 

Ionizing photons from young stars play an important role 
in galaxy evolution and the cosmic reionization of the in- 
tergalactic medium (IGM) in the early universe. Stellar ra- 
diation can ionize interstellar gas, and a fraction of them 
can esca pe from galaxies , and reionize the neutral hydro- 
gen (e.g..llliev et al.ll2007l:lGriedin. Kravtsov. fc Chenll200i 



Razoumov fc Somm cr- L arse nl I2OI0I : lYaiima et al'l 120091 : 
Yaiima. Choi, fc Nagamind 201ll '). The ionization structure 



in IG M is highly relevant to the detectability of high- 2: LAEs 
(e.g., iMcQuinn et all l2007l : iDaval et aLll201ll ). In particu- 
lar, ionization affects significantly the properties of distant 
LAEs, as the emissivity, scattering and escaping of Lya pho- 
tons depend sensitively on the ionization state of the ISM. 

The existing algorithms of ionization RT can be 
broadly divided into three categories: m oments, Monte 
Carlo , and ray-tracing methods (see Trac fc GnedinI 
I2OIII for a rece nt revi ew). Recently, llliev et all ( 20061 ) 
and llliev et all (|2009l ') conducted a comparison of 
various cosmological ionization codes. It was shown 
that the ray-tracing method provides the most accu- 
rate solution but is co m putationally expensive (e.g., 
Gnedin fc Ostrikei 1997: [ Umemura. Nakamoto. fc Susa 



199£ 



Abel. Norman, fc Madau 



199! 



Abel fc Wandelt 



20021 : IRazoumov et all 120021: iMellema et al.l \20oS : ISusa 
2009 : iHasegawa fc UmemurJ bOld). the moments method 
is usuallv diffusive (e.g.. Icnedin fc AbeJ [2003), and the 
Monte Carlo method is efficient and fi exible but requires 



a large number of photon packets (e.g., Ciardi et al.l I2OOII : 



iMaselU. Ferrara. fc Ciardill2003l : [Pierleoni et al.ll2009l ') 

Our implementation of the ionization RT in ART^ code 



is based on the Monte Car lo method, which is similar to 
CRASH (jMaselh et al.ll2003l ). 



2.3.1 Basic Equations and Methods 

The RT of ionizing photons in our code also includes ab- 
sorption and scattering process by hydrogen gas and in- 
terstellar dust. The optical depth is estimated similar to 
equation (fTU)) . dr = nHia-iio{i')ds + mdad(i^)ds , where 
crjjo(!^) = 6.3x W^^^ (u/vl,)'^ is the ioniz ing cross section o f 
hydrogen with Lyman limit frequency ui^ (|Osterbrocklll98S^ ) . 
The time evolution of ionization degree of hydrogen gas is 
estimated by 



dX 
dt 



(21) 



where k^^^^, fcuva photo-ionization efficiency by star 

and UVB, respectively, fc*^ is the coUisional-ionization effi- 
ciency, and «roc is the recombinat ion coefficient. Both k'^ 
and Qrcc are taken from ICenI (| 19921 ). The fc^ by UVB is es- 
timated by self-shielding model, in which UVB is optically 
thin at tih < Wcrit, zero intensity at nu ^ ncrit. We set 
the critical density to be ncrit = 0.0063 cm~^ based on the 
observation d ata of column density distribution of neutral 
hydrogen gas (|Nagamine. Choi, fc Yaumall2010l ). 

The photo-ionization rate by stars and AGNs is es- 
timated by absorbed-photon number in a cell N!^, i.e., 
^star'^Hi = N!^. We treat the time integration with an 
adaptive time step introduced in lBaek et al.l (|2009l ). by us- 
ing three time steps, Atevo, AtaT and Atroc, where Atevo 
is an evolution time step to renew the ionization state. It 
is determined by splitting the recombination time scale of 
volume-mean density by ~ 100, which is close to that of low- 
density region. In each Atevo, the photon- budget of A^'ph are 
traced and the ionization rate is estimated. The evolution of 
ionization degree is renewed with Atevo, which is divided by 
Atrec, a time step estimated from recombination time-scale 
in each cell such that the ionization rate is constant. How- 
ever, when the ratio of the absorbed photon number to that 
of neutral hydrogen in a cell exceeds a pre-set limit {— 0.25 
in this work), the ionization fraction is renewed with At-RT. 

In addition, the frequency of the photon is sampled from 
the source spectrum, similar to that in the continuum RT 
calculation. We use the stellar populati on synthesis model of 
GALAXEV (|Bruzual fc Chariot! I200I) to produce intrinsic 
SEDs of stars for a grid of metallicity an d age, and we use 
a simple , broke n power law for the AGN (|Li et al.ll2008l ). A 
ISalpeteil (|l955l ) initial mass function is used in our calcula- 
tions. 



2.3.2 Test Calculations 

To test our code, we simulate an Hll bubble by a central 
point source in a hydrogen sphere with uniform density, i.e., 
the Stromgren sphere ( jStro mgrcn 1939j). The radius of the 
Hll bubble, r^, is estimated by the balance between the ion- 
ization rate of the source and the recombination rate of the 
gas 



(22) 
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Figure 5. Ionization structure by a point source in a uniform density field. The color bar indicates the neutral fraction of the hydrogen 
gas in log scale. Left panel: dust-free medium with njj = 10""^ cm~^, and right panel: dusty medium with njj = 10~^ cm"'' and 
Tg^j = 4.0, where rg^ is a optical depth of dust absorption from source to the radius of the Stromgrem sphere. The source is placed at 
corner of the simulation box with N\o^ = 5 X 10** and the cell size is 66.7 pc. 



where A^ion is the number of ionizing photons from the source 
per second, and qb is the recombination coefficient to all 
excitation states. 

And the position of the ionizing front of the Hll bubble, 
Ti is analytically derived bv lSpitzen () 19781 ). 

ri = r,{l-e-^y , (23) 

where t^^^ is the recombination timescale. When t 2> free, ^'i 
approaches to rg. 

In the presence of dust, the size o f Hll bubble decreases 
due to dust extinction. ISpitzeil (| 19781 ) has also derived ana- 
lytically the size of the Hll bubble ri in a dusty medium: 

3 y^e^^^-^dy = 1, (24) 
Jo 

where y = r/rs, yi = Ti/rs, and rsd is the optical depth of 
dust from the source to the distance of Stromgren sphere. 

In the tests, we calculate the ionization structure of a 
uniform density field. The source is placed at the corner of 
the simulation box with Mon = 5 x 10** s~^. The number 
density is nn ~ 10^"^ cm~^, and the recombination coeffi- 
cient QB = 2.59 X 10~" cm^ s"^ (T = lO" K). 

Figure [S] shows the neutral fraction of hydrogen gas at 
mid-plane a.tt = 3frec , in a dust-free (left panel) and a dusty 
medium (right panel). The Hll bubble becomes significantly 
smaller than that in dust-free medium. Figure [S] shows the 
corresponding neutral fraction of hydrogen gas as a func- 
tion of the distance from the source for both dust-free and 
dusty tests, and Figure [7] shows the position of the ionization 
front of the dust-free medium. All our res ults are in goo d 
agreements with the analytical solutions of ISpitzeil (|l978l ). 




Figure 6. Neutral fraction of hydrogen gas in log scale as a 
function of distance from source in a dust-free uniform medium. 
The distance is normalized by the radius of Stromgrem sphere. 
The color lines represent our simulation results for a dust-free 
(blue line) and a dusty medium with Tg^ = 4.0 (red l ine), and 
the dotted lines are the analytical solutions bv lSpitzeJ l ll978t) . 



2.4 Adaptive Refinement Grid 

Our RT calculation is done on a grid. In order to handle arbi- 
trary geometry, and cover a large dynamical range while re- 
solving small-scale, high-density gas regions in galaxies and 
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Figure 7. Time evolution of the Hll region radius in a dust- 
fre e uniform med ium. The dotted line is the analytical solution 
by ISpitzeil l ll978l ). and the red filled circles are our simulation 
results. The error bars indicate the radius of neutral fraction of 
0.3 - 0.7. 



ISM, we use an adaptive-mesh refinement grid in 3-D Carte- 
sian coordinate s. The d e tailed description of the adaptive 
grid is given in iLi et al] ()2008h . Here we briefly summarize 
the basic method and procedures. 

We typically start with a base grid of 4"^ box covering 
the entire simulation volume. Each cell is then adaptively 
refined by dividing it into 2^ sub-cells. The refinement is 
stopped if a predefined maximum refinement level, RL, is 
reached, or if the total number of particles in the cell be- 
comes less than a certain threshold (typically set to be 32 for 
SPH simulations), whichever criterion is satisfied first. The 
resolution of the finest level is therefore Lmin = I/box/2^'"'''^, 
where Lbox is the box length, and RL is again the maximum 
refinement level. 

The ART^ code can be easily applied to either grid- or 
particle-based simulations, once the snapshot is converted 
to the format of our grid, which is an effective octo-tree. For 
particle-based simulations, the physical quantities of parti- 
cles are interpolated onto a grid. For example, the gas prop- 
erties at the center of each grid cell, such as density, temper- 
ature, and metallicity, are calculated using the SPH smooth- 
ing kernel of the original simulation. All physical quantities 
are assumed to be uniform across a single cell. 

Figure |8] shows an example of the adaptive grid applied 
to a snapshot redshift z = 3.1 from the cosmological simu- 
lation used in this work (see Section [3. II for details). For the 
parameters used in this application, Lbox = 1 Mpc, RL = 
11, we use a total grid number of only ~ 58300 and reach a 
minimum cell size of Lmin = 59.5 pc in physical scale, which 
corresponds to the spatial resolution of 244 pc in comoving 
scale of the hydrodynamic simulation. This grid resolution 
is equivalent to a ~ 4096'^ uniform grid, which is prohibitive 
with current computational capabilities. 




Figure 8. Example of the adaptive grid applied to a snapshot 
at redshift z = 3.1 from the cosmological simulation on top of 
the gas density distribution of the system. The box size is 1 Mpc 
in comoving scale. The maximum refinement level is 11, and the 
size of the finest cell is 59.5 pc in physical scale. 



2.5 ISM Model 

ART^ implements a m ulti-phase model of the ISM 



(|McKee fc Ostrikeilll977[). as adopt e d in th e hydrodynamic 
simulations ( Springel fc Hernguistl l2003al lbh. It consists 
of a "hot-phase" (diffuse) and a "cold-phase" (dense 
molecular and HI core) components, which co-exist under 
pressure equilibrium but have different mass fractions 
and volume filling factors (i.e., the hot-phase gas is 
^ 0.1 in mass but > 0.99 in volume). The cold gas is 
assumed to follow two empirical relations observed in giant 
molecular clouds: a power- l aw mass distribution (e.g., 
LarsonI 11981): [Fuller fc Mversl Il992l: IWard- Thompson et alj 



19941: [Andre Ward- Thompson, fc MottT ^ Il996l : 

Blitz fc Ro solowsk vll200ir an d a power-law mass-ra dius re- 
lation (e.g. Sanders. Scoville. fc Solomo n 1985; Dame et al 



1986 
2005 



Scoville et al.H 19871 : ISolomon et al.lll987 : iRosolowskvl 



20071 ) . Within each grid cell in the RT calculation, 



the hot gas is uniformly distributed, while the cold, dense 
cores are randomly embedded. 

The dust is assumed to be associated with both the 
cold and hot phase gases, and it follow the distribution of 
the gas with certain dust-to-gas ratio. In the RT calculation, 
when a photon enters a cell, we first determine the distance 
it travels in the hot phase gas before hitting a cold cloud, 
then the distance it travels in the cold cloud before it is 
absorbed. The combination of the two gives the total optical 
depth in the multi-phase ISM, which is then compared with 
a random number to determine whether the photon should 
be stopped for scattering or absorption (see iLi et al.ll2008l 



10 Yajima et al 



for a detailed description and implementation). The Lya 
RT in a multi-phase ISM and dust model is one of the new 
features in our RT code compared with previous works. In 
ART^, depending on the hydrodynamic simulation and the 
nature of the calculation, we can choose either a multi-phase 
or a single-phase model, in which the gas and the dust are 
uniformly distributed with the mean density in a cell, for 
the ISM. 

Galaxies and quasar s at high redshif t s show extinc- 
tion by interstellar dust (jOuchi et al.ll20o3 : iMaiolino et al.l 
1200.1 ) ■ In the early universe, most of the dust is thought 
to be created by Type II supernovae (SNe). Theo- 
retically, several groups have cal culated dust form ation 
in the ejecta of Type II SNe dTodini fc Ferraral 
Nozawa et al. 20031: Schneider. Ferrara. fc Salvaterra 



2001 



Hirashita et al.1 l2005l: iDwck. Gallia no, fc Jones! 120071 '). In 



2004 



particular, iTodini fc Fer rara (2001) have developed a dust 
model based on standard nucleation theory and tested it on 
the well-studied case of SN1987A. They find that SNe with 
masses in the range of 12 — 35 Mq produce about 0.01 of 
the mass in dust per supernova for primordial metallicity, 
and ~ 0.03 for solar metallicity. 

In ART'^, w e ado pt the dust size distribution of 
iTodini fc Ferraral (|200lh for solar metallicity and a M = 
22 Mq SN model, as in Figure 5 in their paper. The size 
distribution is then combi ned with the dust abs orption and 
scattering cross section of I Weingartner fc Drainc. ((2001i ) to 
calculate dust absorption opacity curves. The SN models 
include silicate fractions fs = 0.1, which shows th e silicate 
featur e ~ 9.7 /xm (see the resultant opacity curve iLi et al.l 
l|2008h ). The dust-to-gas mass ratio in each cell is linearly 
proportional to the metallicity in the cell. The ratio is nor- 
malized to Galactic value when the metallicity is the solar 
abundance. 



2.6 Coupling Ionization, Lya Line, and 
Continuum 

To date, there is little work on the couplin g of continuum 
and L ya line transfer. In the seminar work of lPierleoni et all 
(|2009l ), they incorporated pre-computed tables of the prop- 
erties of the Lya photons in the calculation of ionizing UV 
continuum radiation. Our approach differs from theirs in 
that we follow the actual propagation and scattering of the 
LyQ photons in the medium, and take into account the ef- 
fects of ionization and dust absorption. 

In ART^, the photons are treated in three regimes: ion- 
izing continuum at A ^ 912 A, Lya line with center A = 
1216 A, and non-ionizing continuum at 912 A < A ^ lO'^ A. 
We follow the following steps to couple the radiative transfer 
of these three components: 

1. The gas and dust content of each grid cell is determined 
before the RT process. 

2. The neutral hydrogen fraction nni in each cell is deter- 
mined by following the ionizing continuum photons. 

3. The resulting ionization structure is then input to cal- 
culate the RT of Lya photons. 

4. The absorption of the ionizing photons by dust will be 
taken into account in the calculation of the dust re-emission 
in longer wavelengths of the non-ionizing continuum. 



2.7 Monte Carlo Method 

The Monte Carlo RT method follows the propagation, 
scattering, absorption, and reemission of "photon packets" 
(groups of equal-energy, monochromatic photons that travel 
in the same direction), by randomly sampling the various 
probability distribution functions that determine the opti- 
cal depth, scattering angle, and absorption rates. In detail, 
the Monte Carlo ray-tracing procedure involves the follow- 
ing steps: 

1. A photon packet of continuum is emitted from either 
a stellar source or an accreting black hole with random fre- 
quencies consistent with the source spectra. In the case of 
Lya photon, it is emitted from ionized hydrogen gas. The 
photon is emitted with a uniformly distributed random di- 
rection. The probability of a photon being emitted by any 
given source is determined by its luminosity relative to the 
total. 

2. A random optical depth over which the photon must 
travel before an interaction with gas or dust occurs, n = 
— In^, is drawn to determine the interaction location. The in- 
teraction includes scattering and absorption. In our method, 
the photon energies are not weighted, only one event is al- 
lowed. That is, at any given interaction site, the photon is 
either scattered or absorbed, but not both. 

3. Starting from the location of the photon emission, the 
cumulative optical depth of the photon, rtot, is calculated 
stochastically for both hot and cold ISM in a mulit-phase 
model, or single-phase medium with the mean density along 
the ray. If the photon is stopped for interaction within a 
single cell, then rtot is the sum of contributions from possibly 
multiple segments of both hot and cold for the multi-phase 
model, or a single-phase gas/dust within this cell. If the 
photon passes through multiple cells before an interaction 
occurs, then Ttot is the sum of all contributions from relevant 
segments in these cells. 

4. At each boundary between the hot and cold phase gas 
clouds, or at the boundary of the grid cell, the next interac- 
tion point is determined by the comparison between n and 
Ttot- If Ti ^ Ttot, then the photon is either scattered or ab- 
sorbed, with a probability given by the scattering albedo. 
The exact interaction location is then determined inside ei- 
ther hot or cold phase gas, such that rtot becomes exactly 
Ti. If the photon is scattered, its direction is altered using 
a phase function, and the ray tracing of the new photon is 
repeated from step 2. The phase function is the Henyey- 
Greenstein for dust, and simply isotropic for gas. If the 
photon is absorbed by dust, the dust temperature is raised, 
and a new photon is re-emitted according to the scheme de- 
scribed below. The ray tracing of the newly emitted photon 
again restarts from step 2. 

5. If the photon escapes from the system without reaching 
the optical depth Ti, it is then collected in the output spec- 
trum and image. The next photon will be picked up from 
the source, and the whole Monte Carlo procedure from step 
1 will be restarted. 



2.8 Making Images 

In ART^, we incorporate a "peeling off" procedure using 
weighted photons to obtain high-resolution images from a 
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particular viewing direction. From a particular viewing di- 
rection, the contribution from direct (radiation source or 
thermal dust emission) and scattered photons are tracked 
separately with different weight. 

At each location of emission, the contribution from di- 
rect photons to the image is computed as 



code uses the "TreePM" method JXul |l995l') that com- 
bines a "tree" algorithm (|Barnes fc Hut 1986h for short- 



IdiTcc{x,y) = Jo 



47rd2' 



(25) 



where Iq is the energy of each photon packet, x and y are 
the position of the emission site projected onto the image 
plane, and ri and d are the total optical depth and distance 
from the emission site to the observer, respectively. 

At each scattering site, the contribution from scattered 
photons to the image is 



/scatt(a;,i/) = Iq- 



dp 



(26) 



where T2 is the optical depth from the scattering location to 
the observer, and f{6) is the scattering phase function. 

The total intensity of the image is the sum of /direct and 
/scatt, and is accumulated as each photon packet is emitted 
and scattered during the course of the Monte Carlo simu- 
lation. To obtain an image for a given telescope at a given 
waveband, a specific filter function is then applied to con- 
volve the image. 



3 APPLICATION TO COSMOLOGICAL 
SIMULATIONS 

3.1 The Simulations 

The cosmological simulation presented here follows the for- 
mation and evolution of a Milky Way-size galaxy and its 
substructures. The simulation includes dark matter, gasdy- 
namics, star formation, black hole growth, and feedback pro- 
cesses. T he initial condition i s originally from the Aquarius 
Project jSpringel et al.|[2008h . which produced the largest 
ever particle simulation of a Milky Way-sized dark matter 
halo. The hydrodynamical initial condition is reconstructed 
from the original coUionsionless one by splitting each original 
particle into a dark matter and gas particle pair, displaced 
slightly with respect to each other (at fixed center of mass) 
for a regular distribution of the mean particle separation, 
and with a mass ratio correspond ing to a baryon fraction of 
0.16 (|Wadepuhl fc Springej[201ll ). 

The whole simulation falls in a periodic box of 100 
/i~^Mpc on each side with a zoom-in region of a size 
5x5x5 /i~^Mpc^. The spatial resolution is ~ 250 
pc in the zoom-in region. The mass resolution of this zoom- 
in region is ~ 2 x 10^ /i~^Mq for dark matter particles, 
~ 1.9 x 10* H^^Mq for gas and star particles. The cosmo- 
logical parameters used in the simulation are ilm = 0.25, 
Qa = 0.75, (78 = 0.9 and h = 0.73, consistent wit h the 
five-year resuhs of the WMAP (|Komatsu et al.ll2009l '). The 
simulation evolves from z — 127 to z = 0. 

The simulation was performed using the parallel, 
N-body/Smoothed Particle Hydrodynamics (SPH) code 
GADGET -3, which is an improved version of that de- 
scribe d in ISpringel. Yoshida. fc White! (|200ll ) and ISpringell 
(|2005l ). For the computation of gravitational forces, the 



range forces and a Fou rier toansform particle-mesh method 
iHocknev fc Eastwood! [l981i ) for long-range forces. GAD- 
GET implements the entropy-co nserving formulation of 
SPH jSpringel fc Hernguist! !2002!) with adaptive particle 
smoothing, as in !Hernauist fc Katj (|l989! '). Radiative cool- 
ing and heating processes a re calculated assuming coUi- 
sional ionization equilibrium jKatz et al.! !l996! : !Dave et al.! 
!l999h . Star formation is modeled in a multi-phase I SM, with 
a rate that follows th e Schmidt-Kennicutt Law (| Schmidt! 
! 19591 : !Kennicutt|[l99i '). Feedback from supernovae is cap- 
tured through a multi-phase model of t he ISM by an effective 
equati on of state for star-forming gas jSpringel fc Hernguist 
!2003al '). The UV background model of !Haardt fc Madau 
( |l996r ) is used. 

Black hole growth and feedback are also in- 
cluded in our simulation based on the mode l 

of !Springel. Pi Matteo. fc Hernguist! (|2005bh : 

iDi Matteo. Springel. fc HernguistI ^200^ . where black 
holes are represented by coUisionless "sink" particles 
that interact gravitationally with other components and 
accrete gas from their surroundings. The accretion rate 
is estimated from the local g as density and sound speed 
using a spherical Bondi-Hoyle (|Bondi|[l952l : !Bondi fc Hovld 
!l944l : !Hovle fc LvttletonI !l94ll ) model that is limited by 
the Eddington rate. Feedback from black hole accretion 
is modeled as thermal energy, ~ 5% of the radiation, 
injected into s u rround ing ga s isotropically, as descr ibed in 
!Springel et"all l|2005b! ') and !Di Matteo et al.1 (|2005! '). This 
feedback scheme self-regulates the growth of the black 
hole and has been demonstrated to successfully reproduce 
many observed properties of loca l elliptic al galaxies (e.g. 



Springel. Di Matteo. fc Hernguist! !2005a! : !Hopkins et al 



Li et al 



Li et al 



20061) and the most distant quasars at z ~ 6 I 

2007^ . We follow the black hole se eding scheme of 

^2003) and !Di Matteo et al.! (|2008! ) in the simulation: a seed 
black hole of mass Mbh = 10^ h-^MQ was planted in the 
gravitational potential minimum of each new halo identified 
by the friends-of-friends (EOF) group finding algorithms 
with a total mass greater than 10^^h~^ Mq. 

Each snapshot is processed by an on-flying FOE group 
finding algorithm with a dark matter linking length less than 
20% of their mean spacing. Other types of particles are then 
linked to the nearest dark matt er particle. A subs tructure 
detection algorithm SUBFIND jDolag et al.l \200^ ) is then 
applied to each group to search satellite halos. In the sim- 
ulation, each FOF group with its substructure is identified 
as a galaxy (a detailed description of the simulation is pre- 
sented in Zhu et al., in preparation). 

In this work, we apply the ART^ code to the ten most 
massive galaxies in each snapshot at six redshift bins from 
z = 3.1 — 10.2 and focus on the Lya and multi-band proper- 
ties. In our post-processing procedure, we first calculate the 
RT of ionizing photons (A ^ 912 A) and estimate the ioniza- 
tion fraction of the ISM. The resulting ionization structure 
is then used to run the Lya RT to derive the emissivity, 
followed by the calculation of non-ionizing continuum pho- 
tons (A > 912 A) in each cell. Our fiducial run is done with 
Npii — 10^ photon packets for each ionizing, Lya, and non- 
ionizing components. Because the spatial resolution of the 
cosmological simulation is not adequate to resolve the mul- 
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tiple phase of the ISM, we assume a single-phase medium 
in each density grid. The highest refinement of the grid is 
RL = 11. 



3.2 Multi-wavelength Properties of High-redshift 
Galaxies 

3.2.1 Flux Images 

The galaxies in our sample are the 10 most massive ones 
from the cosmological simulation at six snapshots from 
z — 3.1 — 10.2. Note that because our simulation focuses 
on a Milky Way-size galaxy, our sample represents the pro- 
genitors of the Milky Way at different cosmic time, they are 
not the typical "massive" galaxies one expect from a gen- 
eral cosmological simulation with a mean overdensity. The 
total and stellar masses of these galaxies become larger with 
time (decreasing redshift) due to gas accretion and merging 
process. For example, the total mass of the most massive 
galaxy is ~ 2.6 x 10^° Mq, 5.8 x 10^° Mq, 1.6 x 10"Mq, 
and 5.4 x lO" Mq at z = 10.2,8.5,6.2 and 3.1 respec- 
tively, and the corresponding stellar mass is ~ 8.5 x 10* Mq, 

4.3 X 10^ Mo, 9.3 X 10^ M©, and 4.1 x 10^° Mq. 

Figure [9] (left panels) show the projected density of gas 
and stars of each snapshot at the mentioned redshift. Fil- 
amentary structure is clearly seen, and the most massive 
galaxy resides in the intersection of the filaments. The result- 
ing Lya and near infrared (NIR) fluxes of the most massive 
galaxy from our RT calculations are shown in the middle 
and right panels of Figure |9] respectively. The NIR band 
corresponds to IRAC-1 (3.6 fj,m) of the Spitzer telescope. 
The flux image is generated from a "peeling off" technique 
in the RT calculation as described in Section 12.81 

The Lya emission appears to follow the distribution of 
gas and stars of the galaxy. Its surface brightness and size 
increase with time (decreasing redshift), while those of the 
NIR show different pattern. The Lya emission depend on 
the star formation rate (SFR) of the galaxy and the photon 
escape fraction. Although SFR decreases at lower redshift, 
the escape fraction increases as the galaxy grows larger, re- 
sulting in an increase of Lya luminosity. In addition, the 
Lya emission is bright over extended region, owing to many 
scattering processes. 

On the other hand, the observed IRAC-1 in these red- 
shifts corresponds to wavelengths ~ 3200 — 9000 A in the 
rest frame, hence it depends on both SFR and stellar mass of 
the galaxy. The cross section of dust decreases steeply with 
increasing wavelength, it becomes very small at these wave- 
lengths. Therefore, a significant fraction of photons emitted 
by the stars in this wavelength range can escape without 
being scattered or absorbed by dust. Thus, unlike Lya, the 
IRAC-1 image is more compact, and it concentrate on the 
high density peak of stars. We note that the infrared fiuxes 
are higher than those of the LAE sample at z — 3.1 by 
ICronwall et"al] (|200'i1 ). but comp a rable with those of the 
IRAC-detected LAEs in lLai et all lj2008h . This may be due 
to high SFR in these galaxies. 

The NIR wavelength shown here is similar to the 
F356W filter of the next generation telescope, the James 
Webb Space Telescope (JWST). The detection limit of the 
filter in AB magnitude can achieve rn ab = 28.1 with 1 
hour exposure time (e.g. JZackrisson et al.|[201lh . Our calcu- 
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Figure 10. Spectral energy distribution of model galaxies at 
rcdshifts z = 10.2 (green), 8.5 (orange), 7.2 (blue), 6.2 (purple), 
and 3.1 (black), respectively. 



lations show that the model galaxies have mAB = 24.3 {z = 
3.1), 25.1 {z = 6.2), 28.5 (z = 8.5), and 29.8 (z = 10.2) at 
this wavelength. Hence, the JWST will be able to detect the 
infrared emission from galaxies such as the progenitor of the 
Milky Way at z ~ 6 with a ~ 1 minute exposure time. Even 
at higher redshifts, these galaxies may be detected with an 
integration of ~ 2 hours for z — 8.5, and ~ 23 hours for 
z = 10.2. 

3.2.2 Evolution of Spectral Energy Distribution 

The corresponding multi-wavelength SEDs of the most mas- 
sive galaxy at selected redshift are shown in Figure [10] The 
shape of the SED changes significantly from z = 10.2 to 
z = 3.1, as a result of changes in radiation from stars, 
absorption of ionizing photons by gas and dust, and re- 
emission by dust in the infrared. The Lya line appears 
to be strong in all cases. The deep decline of UV contin- 
uum at 2 > 8 is caused by strong absorption of ionizing 
photons by the dense gas. Galaxy at lower redshift has a 
higher floor of continuum emission from stars and accret- 
ing BHs, a higher ionization fraction of the gas, and a 
higher infrared bump owing to increasing amount of dust 
and absorption. Moreover, due to the effect of negative k- 
correction, the flux at > 500 ^m stays close in different 
redshifts. Our calculations show that the model galaxies 
have a flux of = 0.043 (2 = 3.1), 0.057 (2 = 6.2), 
0.02 (2 = 8.5), and 0.004 (2 = 10.2) mjy at 850 fj.m in ob- 
served frame. The new radio telescope, Atacama Large Mil- 
limeter /submillimeter Array (ALMA) may be able to detect 
such galaxies at 2 ~ 6 with ~ 2 hours integration, and at 
2 ~ 8.5 for ~ 20 hours with 16 antennas. Even with ALMA, 
it may be very difficult to detect such a galaxy at 2 = 10.2 
(more than 70 days integration). 

From the SEDs, it is interesting to compare the emer- 
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Figure 9. Evolution of galaxies in different wavebands. The left panels show the distribution of gas and stars from the cosmological 
simulation. Middle panels show the surface brightness of Lya, and right panels are the flux images in near infrared IRAC-1 (3.6 iim) on 
board of Spitzer telescope. The box size is 1 Mpc in comoving scale. Note the middle and right panels show only images of the most 
massive galaxies corresponding to the rodshifts labeled in the left column. 
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Figure 11. Lyo properties of modeled galaxies at different red- 
shifts, in clock wise direction: star formation rate, emergent Lya 
luminosity, equivalent width of Lyo line in rest frame, and photon 
escape fraction of Lyo (open circles) and UV continuum (1300 - 
1600 A, crosses). Red and Blue filled circles are the values of 
the most massive galaxies and the median in our sample at each 
snapshot, respectively. Filled and open red circles in the panel of 
escape fraction shows the Lya and UV continuum escape fraction 
of the most massive galaxies, filled and open blue circles are the 
median values of Lya and UV continuum, respectively. 



gent Lya luminosity and that of the infrared integrated from 
8 to 1000 /xm in the rest frame. A simple relation between in- 
frared and Lya luminosity has been suggested bv lKennicuttI 
l|l998l) through star formation, namely SFR(Mq yr~^) — 
4.5 X IQ-** LiR(ergs s"^) and SFR(Mq yr"^) = 0.91 x 
lO"*^ Lc«(ergs s~^) with the assumption of La/Lua = 8.7 
(cEise B). However, our calculation suggests a more complex 
relation between Lya and infrared luminosity, as inferred 
from Figure [TOl This is because La depends on ftZ", while 
the LjB. depends on f^^Y, and the two escape fractions differ 
from each other. In our simulations, the fHY mostly de- 
creases with redshift, and is below 0.5 at 2 < 6. Hence, the 
UV continuum may serve as a good tracer of SFR at z > 6, 
while 2 < 6 the infrared flux may be one. In addition, since 
the modeled galaxies at z < 6 show a smaller value of /^c" 
than 0.5 (see next subsection), hence Lya flux also may not 
be a good indicator of SFR. 



3.2.3 The Lya PropeHies 

The derived properties of Lya emission from these modeled 
galaxies are shown in Figure llll in comparison with their 
corresponding SFR. The most massive galaxies at very high 
redshift {z > 6) maintain a high SFR of ~ 10 - 30 Mq yr"^ 
fueled by accretion of cold gas and merging process. The 
SFR decreases to ~ 3 M0 yr~^ at z ~ 3 owing to depletion 



of star-forming gas and suppression by radiation feedback 
from both stars and accreting BHs. The resulting emergent 
Lya luminosity from recombination of ionizing photons and 
from excitation cooling increases from ~ 1.0 x 10*^ ergs s~^ 
at 2 = 8.5 to ~ 2.3 x 10*^ ergs s"^ at 2 = 3.1. The observed 
LAE luminosity functions show a chara cteristic luminosit y 
4.4 X 10"^ ergs s"^ at 2 = 6.6 JOuchi et aLlbOldl . 
6 X 10^ ergs s~ ^ at 2 = 3.1 jOuchi et all ' 



20081 : 



CSC of Lya and 
~ 0.75 



iGronwall et al.ll2007l : [b1 anc et al]l201ll ). which are compa- 
rable to the luminosity range of our modeled galax;ies. The 
median values of SFR moderately decrease with increasing 
redshift, and show ~ 0.5 — 3 M0 yr~^. However, the La 
does not decreases largely by the redshift evolution of /esc. 
In addition, a few galaxies in each snap shot are Lya bright, 
and have Lya luminosity oi La 10^^ ergs s~^, which is de- 
tecta ble with recent observational studies (e.g.. lOuchi et all 
l2008i ). 

The lower-left panel of figure [TT] shows the photon es- 
cape fraction of Lya, /^c" a-nd UV continuum /^J, where 
ffsY is calculated at Arcst = 1300 - 1600 A. Both escape 
fractions appear to increase with time (decreasing redshift), 
with /c^sc" ~ 0.18 at 2 = 8.5 to ~ 0.5 at 2 = 3.1, while 
ffj ~ 0.09 at 2 = 8.5 to ~ 0.6 at 2 = 3.1. In the early 
evolution phase, the most massive galaxies are compact and 
gas-rich. They grow rapidly from accretion of cold gas which 
fueled intense star formation. The highly concentrated gas 
and dust efficiently absorb the Lya and UV photons, and 
suppress their escape. At a later time, the gas and stars dis- 
tribute to a more extended region, and a large fraction of 
gas is converted to stars. Hence, the photon may easily es- 
cape from the galaxy, contributing to a higher photon escape 
fraction. 

On the other hand, the median value of 
UV increases monotonically with redshift, e.g., /esc 
at 2 = 8.5 to ~ 0.29 at 2 = 3.1, while ffj ~ 0.81 at 
2 = 8.5 to ~ 0.3 8 at 2 = 3.1. Such a trend is consistent with 
observations by iHaves et all (|201ll ). 

Unlike the most massive galaxies, the low-mass ones 
have lower SFR and lower dust content, the ISM is less dense 
and less clumpy, and the metallicity is lower. As a result, the 
escape fractions are higher in small galaxies than the massive 
ones. 

The EW of Lya line is defined by the ratio between the 
Lya fiux and the UV fiux density /uv in rest frame, where 
the mean flux density of A = 1300 — 1600 A in rest frame is 
used. In practice, the EW is estimated by EW — EWint ffv ; 

GSC 

where EWint is the intrinsic equivalent width. The EWint 
depends on the stellar SED and the excitation Lya cooling: 
the younger SED (or top-heavy IMF) and the more efficient 
Lya cooling, the higher EWint. 

The resulting Lya EWs of the most massive galaxies are 
shown in the lower right panel of Figure lTT] Since the photon 
escape fraction of Lya is close to that of the UV continuum, 
the EW shown is basically EWint. All three modeled galax- 
ies h ave EW > 20 A, they are therefore classified as LAEs 
(e.g.. IGronwall et al.ll2007l ). The EW evolves gradually with 
redshift. At 2 > 6, the EWs are high, EW ~ 60 A, owing 
to strong Lya emission from excitation cooling and AGN 
activity, while at 2 = 3.1, it drops to ~ 21 A as the Lya 
emission is mainly produced by stellar radiation. The me- 
dian value increases more steeply with redshift, from ~ 14 A 
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Figure 12. Emergent Lya profile. Dot and soUd lines are the 
profiles of intrinsically emitted and escaped photons respectively. 



at z = 3.1 to ~ 164 a± z — 8.5. These EWs are within the 
observed range, and the trend is in broad agreement with 
observations that galaxies at higher redshifts appear to have 
higher EW than their counterparts at l ower redshifts (e.g., 
iGronwall etHI l2007l : lOuchi et ahllioOSl ). Some galaxies at 
higher redshift z > 8 show very high EW because the exc i- 
tation Lya cooling becomes dominant (|Yauma et al.|[2011bl ). 
However, it will be suppressed significantly by IGM at such 
a high redshift. 

The recent discove ry of the most distan t LAE, UDFy- 
38135539, aX z = 8.6 bv lLehnert et al.1 (|2010t ) marks a mile- 
stone in observational cosmology that galaxies form less than 
600 million years after the Big Bang. The detected Lya line 
flux of this object is ~ 6 x lO"'^** ergscm"^ s~^, which im- 
plies a Lya luminosity of ~ 5.5 x lO*'^ ergss"^. This is close 
to that of our galaxy ai z = 8.5. Our model suggests that 
this LAE may have a total mass of ~ 10^"^^^ Mq, a stellar 
mass of ~ 10** Mq, a SFR of the order of 10 M p yr'^ and 
fcZ'^ ~ 0.18. We note that iDaval et al.1 (120111 ) suggested a 
lower SFR of 2.7 — 3.7 M© yr~^ for this galaxy by assum- 
ing a Lya escape fraction of 100 per cent. However, our de- 
tailed RT calculations show that the escape fraction is much 
smaller than unity, which means it would require a higher 
SFR in order to produce the observed Lya flux. 



3.2.4 The Lya Line Profiles 

Figure [T^] shows the profiles of Lya line of the modeled 
galaxies. The frequency of the intrinsic Lya photon is sam- 
pled from a Maxwellian distribution with the gas tempera- 
ture at the emission location in the rest flame of the gas. We 
collect all escaped Lya photons and estimate the line proflle. 
In practice, the inhomogeneous gas structure can cause the 
change of flux depending on the viewi ng angle. Although th e 
Lya flux can change by some factors (lYaiima et al.ll2011bD . 
we focus on only the mean value in this paper. 



The Lya lines of our galaxies show mostly a single peak, 
a com mon profile of LAEs at high redshifts (e.g.. lOuchi et all 
l2010t ). In a static and optically thick medium, the profile 
can be double peaked, but when the effective optical depth 
is small due to high relative gas speed or ionization state, 
there might be only a single peak (|Zheng fc Miralda-Escudg 
l2002f ). In our case, the fiow speed of gas is up to ~ 
300 km s~^, and the gas is highly ionized by stellar and 
AGN radiation, which result in a single peak. 

In addition, the profile at higher redshift shifts to 
shorter wave length, and it shows the chara cteristic shape 
of gas inflow (|Zheng fc Miralda-Escudil2002l ). Ahhough our 
si mulation includes feedb ack of stellar wind similar to that 
of ISpringel et al.l (|2005bh . the Lya line profile indicates gas 
infiow or symetric in these galaxies. Our result suggests that 
high-redshift star-forming galaxies may be fueled by efficient 
infiow of cold gas from the filaments. We will address this 
issue in detail in a forthcoming paper (Yajima et al. in prepa- 
ration) . 

We note that there is a number of studies on inflow and 
outflow in Lya lines, both obser vationally (iKunth et al.l 



1998h and theore t ically ( e.g., Verhamme et al.l 



Diikstra fc Wvithe 



2006 



2OIOI : iDiikstra. Mesinger. fc Wvithd 



20 111 : iBarnes et alllionT It has been suggested that the 



IGM correct ion of Lya emission may be en hanced due to 
inflow (e.g., IDiikstra. Lidz. fc WvitQ |2007^ or decreased 
via outflows fe.g.. IVerhamme et al.ll2006t ). However, there is 
a large uncertainty in the prescription of wind feedback in 
current simulations. To investigate the velocity fleld of ISM 
in high-redshift galaxies, one needs simulations with a real- 
istic wind model and a high spatial resolution comparable 
to that of line observations, which are currently unavailable 
and are beyond the scope of this paper. Moreover, we 
should point out that the line proflle can change due to the 
difference of the sub-grid ISM models. 

In addition, the line proflle may be highly suppressed 
and change d by scatter i ng in the intergalactic medium 
(IGM ) (e.g., ISantosll2004l: IDiikstra et al.|[2007l: IZhetig et"al] 
I2OIOI : iLaursen. Sommer-Larsen. fc Razoumovl 201 ll ). The 
IGM effectively absorbs the Lya photons at the line cen- 
ter and at shorter wavelengths by the Hubble flow (e.g., 
iLaursen et al.l 1201 ll ). Therefore, the inflow feature in our 
proflles may be disappeared and the shape may change to an 
asymmetric single peak with only photons at red wing. We 
will investigate the IGM correction in future work that in- 
cludes propagation and scattering of Lya and ionizing pho- 
tons in the IGM. 



3.2.5 The Properties of Ionizing Photons 



Ion 

CSC 1 



Figure fTiD shows the escape fraction of ionizing photons, /, 
and the resulting emissivity from both stars and AGNs of the 
modeled galaxies. The escape fraction of the most massive 
galaxies has a range of fl°^ ~ 0.03 — 0.37, consistent wit h 
the resuhs of IDiikstra et al.l (120111 ): lYaiima et al] (|2011al ). 
In addition, fl°^ decreases with increasing redshift. Unlike 
the Lya and non-ionizing UV photons, the ionizing photons 
can be absorbed by hydrogen gas and dust. At high red- 
shift, the star formation efficiency is small, so a large frac- 
tion of hydrogen gas remains neutral and prevent the escape 
of ionizing photons. The median values do not decrease with 
redshift due to lower dust mass. 
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Figure 13. The escape fraction of ionizing ptiotons (top panel) 
and the resulting emissivity (bottom panel) of modeled galaxies 
at redshifts z = 3.1 — 10.2. Red and Blue filled circles are the 
values of the most massive galaxies and the median in our sample 
at each snapshot, respectively. 



Figure 14. Effects of AGN on escape fraction of Lyo photons 
(top panel) and the emergent luminosity (bottom panel). The ion- 
ization and Lycf luminosity are estimated with radiation from the 
following sources: UVB only (triangles), UVB + star (squares), 
and UVB + AGN + star (circles). 



The emissivity of the ionizing photons roughly in- 
crease with decreasing redshift, reaching ~ 5 x 10^^ at 
z = 3.1. in this work, we do not calculate the RT in the 
IGM, but the size of the Hll bubble around the galaxy 
may be estimated by a modified Strom grem sphere. As- 
suming a dumpiness factor of C = 5 (|lliev et al.l l2007l : 
IPawlik. Schave. fc van Scherpenzeelll2009l ), the radius of the 
Hll bubble is estimated to be r = 554, 198, 47 and 28 kpc 
at 2 = 3.1, 6.2, 8.5 and 10.2, respectively. The Lya photons 
may pass through the Hll bubbles, resulting in a correlation 
be tween high-redshift LA Es and Hll regions, as suggested 
bv lMcQuinn et all (|2007l ). In addition, the median value of 
ionizing-photons emissivity is smaller due to lower SFR, and 
hence the Hll bubble becomes smaller than that of the most 
massive galaxy by a factor of ~ 2 — 6. 



4 DISCUSSION 

In observations of Lya emission from high-redshift galaxies, 
there are a number of uncertainties rega r ding the contribu- 
tions from AGNs (e.g. IWang et al]|2004l : lOuchi et al1l2008l : 
iNilsson fc M0lleill2Olll '). stars, and excitation cooling from 
cold accretion. With our model, we can separate the ra- 
diation components and disentangle the contribution from 
different sources. Also, in numerical simulations, resolution 
plays an important role in the robustness of the results. So 
we discuss these physical and numerical effects on Lya emis- 
sion in this Section. 



4.1 Effects of AGN and Stellar Radiation on Lya 
Emission 

The effects of AGN and stars on the Lya emission is shown 
in Figure [21 It appears that the presence of AGN does not 
change significantly the Lya escape fraction. It is because 
AGNs are usually located in galactic center where the gas 
density peaks, most of the ionizing photo ns near Lyman- 
limit 13.6 eV) are tra pped there (e.g., iNagamine et al.l 
l20ld : lYaiima et al.|[2oTlal ) and cannot escape. On the other 
hand, a fraction of high energy photons can escape from the 
galaxy without ionizing interstellar gas due to small cross 
section of gas and dust. Therefore, AGN radiation has little 
effect on the /esc of Lya photons. There is no change in /esc 
in the case where there is no stellar and AGN radiation. 

In addition, the AGN does not affect the Lya luminos- 
ity, as shown in the bottom panel of Figure [TJ] This may be 
because the BHs in these galaxies are still quite small and 
the accretion rates are quite low. On the other hand, the 
presence of stars appears to change the escape fraction and 
emergent luminosity of Lya, because the stellar radiation 
is much stronger than that of AGN in the model galaxies, 
and the effect becomes stronger at lower redshift. Because 
Lya emissivity is produced by the recombination of ioniz- 
ing photons and coUisional excitation from hydrogen gas, 
the interplay between these two mechanisms determines the 
outcome. At redshifts z > 6, excitation cooling from cold 
accretion of dense gas becomes comparable to recombina- 
tion cooling, while at 2 ~ 3, the number of Lya photons by 
recombination in Hll regions by stars and AGNs outnum- 
bers that from excitation cooling. Therefore, the contribu- 
tion from stars to the Lya emission changes with cosmic 
time and environment. The Lya emissivity from excitation 
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Figure 15. Effects of plioton number on the escape fraction of 
hya piiotons. Open triangles, squares and circles indicate photon 
packet number Afp^ = 10^, 10'^, 10*, and 10^, respectively. The 
fiducial number used in our general RT calculations is A'^ph = 10^. 



coolinK can be reduced due to the stellar ionization (e.g., 
iFaucher-Giguere et all 1201(1 ). Moreover, the distribution of 
excitation Lya cooling can be more extended than star- 
forming regions, resulting in higher escape fraction. Thus, 
the escape fraction and emergent luminosity of Lya some- 
what increase in the case without stellar radiation. 

4.2 Resolution Effect 

Figure [T^ demonstrates effects of photon number on the es- 
cape fraction of Lya photons, with different number of pho- 
ton packets, A'^ph = 10^, 10'^, 10* and 10^. It appears that the 
results converge once the photon number is above 10*, which 
is an order of magnitude lower than the fiducial A'ph ~ 10^ 
in our general RT calculations. The relative difference be- 
tween A'ph ~ 10* and 10^ is 1-6 per cent. Of course, the 
photon number should increase with galaxy mass, but for 
our modeled galaxies, A'ph — 10^ should be sufficient. 

The effects of grid resolution on the escape fraction of 
LyQ photons is presented in Figure 1161 where results from 
different maximum refinement level (RL) are compared. It 
appears that /esc changes dramatically with RL, but it con- 
verges at RL J5 10. Grids with poor resolution can produce 
artificially high /esc by a factor of more than 20. Because 
most of the Lya photons are created in high-density gas 
clumps, where they are also effectively absorbed by numer- 
ous scatterings, sufficiently resolving the gas density field is 
important in the calculation of Lya RT. If we follow the 
resolution of RL = 10 with uniform grids, it would require 
a number of cells of ~ 2048'^, too expensive for RT simu- 
lations with the current computational facility. In fact, the 
current limit is ~ 500'^, as demonstrated by recent work on 
three-dimensional RT calculations using uniform grid (e.g.. 
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Figure 16. Effect of grid resolution on the escape fraction of Lya 
photons. Different color indicates different maximum refinement 
level (RL): purple - RL = 4, blue - RL = 6, green — RL = 8, 
orange - RL = 10, and red - RL = 11. The fiducial RL used 
in our general RT calculations is 11. The cell size of the highest 
refinement level is 1/2^^^^ Mpc in comoving scale. 

Illiev et al.ll2007l : lYaiima et all [20 Hal ). Therefore, adaptive 
refinement grid is critical in the RT calculation of Lya emis- 
sion. 



5 SUMMARY 

We have presented an improved, three-dimensional, Monte 
Carlo radiative transfer code ART^ to study multi- 
wavelength properties of galaxies. It has the following es- 
sential features: 

• It couples Lya line, ionization of neutral hydrogen, and 
multi-wavelength continuum radiative transfer. The Lya 
module includes emission from both recombination and colli- 
sional excitation, and the continuum includes emission from 
X-ray to radio. 

• It employs an adaptive grid scheme in 3-D Cartesian 
coordinates, which handles an arbitrary geometry and cov- 
ers a large dynamic range over several orders of magnitude. 
Such a grid is critical in resolving the clumpy, highly inho- 
mogeneous density field in the RT calculations. 

• It adopts a two-phase ISM model in which the cold, 
dense clouds are embedded in a hot, diffuse medium in pres- 
sure equilibrium. In the case where hydrodynamic simula- 
tions have sufficient resolution to resolve the multi-phase 
ISM, this model ensures an appropriate prescription of the 
ISM physics, which is important for studying dust properties 
in galaxies. 

• It includes a supernova-origin dust model in which the 
dust is produced by Type-II supernovae, which is especially 
relevant for dust in high-redshift, young objects, because 
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there are insufficient low-mass, old AGB stars in these sys- 
tems to produce the observed amount of dust from the clas- 
sical dust models. 

• It follows the radiation from stars and AGNs, and cal- 
culates the scattering, absorption, and re-emission of the 
surrounding medium. It efficiently and self-consistently pro- 
duces SEDs and images in a full spectrum, including Lya 
emission line and continuum from X-ray to millimeter, as 
well as ionization structures of the medium, for direct com- 
parison with multi-band observations. 

• It can be easily applied to either grid- or particle-based 
hydrodynamical simulations and has broad applications in 
cosmology. 

We have tested the new ART^ extensively, and found 
that the coupling of Lya line and continuum enables a self- 
consistent and accurate calculation of the Lya properties of 
gala^xies, as the equivalent width of the Lya line depends on 
the UV continuum, and the escape fraction of Lya photons 
strongly depends on the ionization structure and the dust 
content of the object. 

We applied the code to a cosmological SPH simulation 
of a Milky Way-like galaxy, and studied in details a sam- 
ple of massive galaxies at redshift z=3.1 - 10.2. We find that 
most of these galaxies are Lya emitters. The escape fraction 
of Lya photons of the most massive galaxy in each snap shot 
increeises with decreasing redshift, from about 0.18 at red- 
shift z ~ 8.5 to 0.6 at z ~ 3. The emergent Lya luminosity 
shows similar trend, increasing from La — 1.0 x 10''^ ergs s"^ 
at 2 ~ 8.5 to La = 2.3 x 10*^ ergs s~^ at z ~ 3. The equiv- 
alent widths of these galaxies change with redshift as well, 
from ~ 67 A at 2 ~ 8.5 to ~ 21 A at 2 ~ 3. The pro- 
files of the resulting Lya lines commonly show single peak, 
and the peak of the profile at 2 > 6 is shifted to shorter 
wavelength, implying gas infiow from surrounding filaments. 
On the other hand, the median values in our gala^cy sample 
shows somewhat different trend from that of the most mas- 
sive ones. The escape fraction monotonically increases with 
redshift, while the Lya luminosity shows a weak evolution, 
as a result of the interplay between SFR and escape fraction. 

Our results suggest that the LAEs and their properties 
evolve with cosmic time, and that the first LAEs in the early 
universe, such as currently th e most distant one detected at 
2 ~ 8.6 (jLehnert et al.|[2010l ). may be dwarf galaxies with 
a total mass of ~ 10^""^^ Mq fueled by cold gas accretion. 
They may have a star formation rate of order of 10 M0 yr~^, 
and a Lya escape fraction of ~ 20 per cent. 
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